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Abstract. - We investigate flow dynamics in rivers characterized by basin areas and daily mean 
discharge spanning different orders of magnitude. We show that the delayed increments evaluated 
at time scales ranging from days to months can be opportunely rescaled to the same non-Gaussian 
probability density function. Such a scaling breaks up above a certain critical horizon, where a 
behavior typical of thermodynamic systems at the critical point emerges. We finally show that 
both the scaling behavior and the break up of the scaling are universal features of river flow 
dynamics. 
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<— I I One of the grand challenges in hydrology is understand- 
1 i ng the real nature of river flow dynamics. In fact, fluctua- 
tions in river flow are the result of complicated interactions 
among several external factors, such as climatic changes 
and landscape characteristics. Several attempts have been 
OOmade to capture either stochastic or deterministic signa- 
t'ures in the dynamics of fluctuations at different time and 
; space scales. However, a definitive comprehension of the 
phenomenon is still missing. 
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c^) \ River flow clearly shows long-range correlations and 
anomalous scaling typical of multifractal processes [IJ[2] ■ 
The Bramwell-Holdsworth-Pinton (BHP) distribution, de- 
k> scr ibing fluctuations of global quantities in correlated tur- 
rNljulent and critical phenomena [31 H], has been recently 
^proposed to explain universal scaling signatures for fluc- 
tuations in river level (5] and river flow [6], after rescaling 
time series with respect to the seasonal mean and variance. 
However, the nature of the fluctuations and multifractal 
behavior can not be uniquely related to a stochastic or a 
chaotic deterministic dynamics, because the observed dy- 
namics changes from a less complex (deterministic) to a 
more complex (stochastic) one by varying the scale of ag- 
gregation [THTTj. Thus, apparently controversial results 
emerge from the analyses of river flow databases, because 
of the observed transitions from determinism to stochas- 
ticity with increasing time scale (see [Hj and Ref. therein) . 



systems. The key point of the analysis is to consider the 
time series of increments, a set of N new time series, each 
constructed as the difference between the original time 
series, locally detrended, and a copy of it delayed by n 
(n = 1,2,...,N) time units. By investigating the proba- 
bility density functions (PDFs) of the increments, we find 
a scaling of the non-Gaussian curves obtained at different 
time scales, that is typical of processes with anomalous 
diffusion. Moreover, we reveal the universality of such 
a scaling by comparing several different rivers. Our re- 
sults suggest the presence of turbulence-like behavior in 
the fluctuations of river flows at any time scale, from daily 
to monthly, up to a certain time horizon, above which a 
behavior typical of thermodynamic systems at the critical 
point emerges. 

We analyze historical time series of the average daily 
discharge as mesaured at several sites in the northeast- 
ern gulf of Mexico (NEGOM) region (USA), that covers a 
large area corresponding to a longitude ranging from —91° 
to -82° and a latitude ranging from 28° to 33°. More- 
over, we consider several other rivers around the world, 
from Asia, Europe and Australia. The list of all sites used 
for this study is shown in Tab.Q] time scries of different 
size, measured in different epochs, and corresponding to 
stations with different physical features, as territory and 
drainage area, are present. A representative example of 
the apparently irregular behavior of river flow, is given 
in Fig.[IJa), where ten years of measurements of the Al- 
abama river discharge, at the Millers Ferry station (AL), 



In this Letter we present an analysis of river flow fluc- 
tuations based on a method originally introduced for the 
investigation of turbulent flows and complex hierarchical 
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are shown. In Fig.QTb) we report the series of delayed and 
detrended increments for the scale n = 2 days, obtained 
as successively explained in the text. 

We firstly investigate the average daily flow q K , k = 
1, 2, 34, of all stations in Tab.[TJ By labelling with q K (t) 
the measurement at time t, t = 1,2, ...,N K days, at the 
k — th station, we inspect both the autocorrelation function 
and the power spectrum of each data set. The former is 
defined from the autocovariancc 

N k -t 

C«(r) = [(«*(*) - <?*)) + T ) - 

K t=i 

as C K (r)/C K (0), and quantifies the amount of correlation 
among measurements delayed by an interval of amplitude 
r (the variance of the series corresponds to C K (0)). The 
delay r* such that the autocorrelation reaches zero for 
the first time, is known as decorrelation time, and esti- 
mates the maximum range of dependence among samples. 
The power spectrum can be defined through the Wiener- 
Khinchin theorem as P K {v) = J r [C K (r)/C K (0)], where T 
is the Fourier transform, which quantifies the amount of 
correlation among samples in the frequency domain. In or- 
der to reduce fluctuations at higher frequencies, the power 
spectrum of a single data set is obtained by dividing q K 
into 16 non-overlapping subsamples, by calculating the pe- 
riodogram for each subsample, and by estimating the en- 
semble average. As a representative example, in Fig.^a) 
is shown the autocorrelation function for all rivers corre- 
sponding to k = 1,2,..., 17. Time series have been pre- 
viously deseasonalized by rcscaling with respect to the 
seasonal mean and variance, as in Refs. [B1[T31[T1]. Al- 
though the considered rivers are located in the same re- 
gion (namely, NEGOM) , the autocorrelation analysis does 
not reveal features common to the different series ana- 
lyzed. In particular, we observe slowly decay correlations, 
with different decorrelation times t* ranging from 70 days 
to 250 days. In Fig.|5fb) we show the logarithm of the 
power spectrum vs. the logarithm of the frequency for 
the same rivers: linear regions, corresponding to scaling 
regimes P K (y) ~ v" k are evident, and can be characterized 
by the scaling exponents f3 K . Scaling regions of different 
length are present on different frequency ranges. In or- 
der to perform a simple comparison among the considered 
rivers, we consider only three particular frequency ranges 
and, for each range separately, we estimate the scaling ex- 
ponents (3 K for each river by considering linear regions of 
the same length. In Fig. [51(c) , we show the rank-ordered 
scaling exponents /3 K , extracted from the slopes of the 
linear regions in Fig.|2Jb), for log(i//day _1 ) 6 [—2,-1], 
[—3,-2] and [—4,-3]. The main result in the frequency 
domain is that, in general, time series measurements for 
each river flow are correlated, although they are charac- 
terized by several scaling regimes with exponents ranging 
from about -0.5 to about -5.0, and corresponding, in prin- 
ciple, to dynamical systems of very different nature. 
A more advanced approach to the study of correlations 



in time series involves the detrended fluctuation analy- 
sis (DFA) [T5] or its generalized version, the multifractal 
DFA (MDFA) [Tg]. The MDFA is a powerful tool, strictly 
related to standard multifractal analysis, allowing i) to 
explore the scaling relationship between detrended fluc- 
tuations and time as a function of a real variable 7, and 
ii) to characterize correlations by scaling exponents that 
generalize the Hurst parameter, the quantity widely used 
to quantify the amount of persistence or anti-persistence 
in time series. The Hurst parameter T-L K , obtained from 
MDFA when 7 = 2, is strictly related to scaling exponents 
(3 K , obtained from spectral analysis, by — /3 K = 2TL K — 1 

An accurate estimation of the Hurst parameter, through 
both MDFA and wavelet techniques, for several rivers 
around the world has been recently reported [T3l[T4l[T5] . 
The analysis put in evidence that detrended fluctuations 
are characterized by different values of the Hurst param- 
eter. In particular, % varies from river to river, ranging 
from 0.55 to 0.95, underlying the long-term persistence of 
river flow dynamics on time scales from weeks to decades 
but no universal behaviour |13j . Moreover, the multifrac- 
tal analysis revealed a self-afflne scaling behaviour and a 
universal function, derived from a generalization of the 
multiplicative random cascade model, has been found to 
describe the multifractal spectra of all rivers with extraor- 
dinary precision [T3J[TH]. River flow fluctuations should 
be strongly dependent on rainfall dynamics, and therefore 
multifractal cascade models should be the best candidates 
for their characterization [T]. However, it has been re- 
cently shown that precipitations might be not responsible 
for both the long-term correlations observed in river flow 
data and the appropriateness of multiplicative cascades 
for their modelling [TJ]. Although such results represent 
an important advance in understanding the dynamics of 
river flow fluctuations, a universal scaling behaviour in the 
time domain is still missing. 

The goal of our spectral study is to show the great vari- 
ance in the scaling exponents f3 K estimated from three dif- 
ferent frequency intervals and corresponding to time scales 
ranging from few days to few weeks. In order to compare 
our results with those reported in literature, for instance 
in Tab. 1 of Ref. [T3], we estimate the power spectrum of 
deseasonalized time series and the corresponding values of 
fi K for —7.6 < log(zy/day -1 ) < —4.5 in the frequency do- 
main or, equivalently, for 90 < n/day < 2000 in the time 
domain. We have chosen this interval because it is compat- 
ible with the interval where Hurst parameters, reported in 
Ref. [13) . have been determined for several rivers around 
the world. For instance, the time series for Thames river 
at Kingston (k = 20) is common to both analyses: in this 
case we obtain /3 2 o = —0.70, corresponding to the value 
W20 = 0.85 ± 0.08 for the Hurst parameter, in good agree- 
ment with the value h{2) = 0.8 reported in literature 
However, in this Letter, we are not interested in a direct 
comparison between our results and those ones obtained 
from MDFA |T31[TB] , but we are interested in investigat- 
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Station name 


IVIean discharge 


Historical data 






(m /aay) 


(irom - 


to) 


1 


Alabama (Millers Ferry, AL) 


80300300 


10/1969 


- 09/2000 


2 


Apalachicola (Sumatra, FL) 


63265600 


10/1977 


- 04/2001 


3 


Fombigbee (Demopolis, AL) 


57410700 


10/1927 


- 04/2001 


4 


Pascagoula (Graham Ferry, MS) 


25518700 


10/1993 


- 04/2001 


5 


Aucilla (Scanlon, FL) 


17600000 


08/1976 


- 09/1998 


6 


Choctawhatchee (Bruce, FL) 


17101400 


10/1930 


- 04/2001 


7 


Suwannee (Branford, FL) 


16804200 


10/1930 


- 04/2001 


8 


Pearl (Monticello, MS) 


16548000 


10/1938 


- 04/2001 


9 


Escambia (Century, FL) 


15348000 


10/1934 


- 04/2001 


10 


Escambia (Molino, FL) 


14333400 


10/1982 


- 03/2000 


11 


Perdido (Barrineau Park, FL) 


1903000 


10/1940 


- 09/2000 


12 


Ochlockonee (Bloxham, FL) 


4148310 


10/1925 


- 04/2001 


13 


Yellow (Milligan, FL) 


2742960 


10/1937 


- 04/2001 


14 


Blackwater (Baker, FL) 


832902 


10/1949 


- 04/2001 


15 


Steinhatchee (Cross City, FL) 


763213 


10/1949 


- 01/2000 


16 


Econfina (Perry, FL) 


342289 


10/1949 


- 09/2000 


17 


Fenholloway (Foley, FL) 


124025 


10/1955 


- 09/1999 


18 


Fish Creek (Ketchikan, AK) 


1037462 


12/1938 


- 12/2009 


19 


Inabon (Real Abajo, Puerto Rico) 


45771 


07/1964 


- 12/2009 


20 


Fhames (Kingston, UK) 


5634413 


01/1883 


- 12/2009 


21 


Severn (Bewdley, UK) 


5252122 


04/1921 


- 12/2009 


22 


Fana (Polmak, Norway) 


14535168 


02/1946 


- 12/1999 


23 


Namsen (Bjornstad, Norway) 


3266237 


10/1934 


- 12/1997 


24 


Verdalsvassdraget (Grunnfoss, Norway) 


3297774 


11/1951 


- 12/2000 


25 


Chiang Mai (Ping River, Thailand) 


5016590 


01/1921 


- 12/1999 


26 


Nakhon Sawan (Chao Phraya, Thailand) 


60307717 


01/1906 


- 12/1999 


it i 
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7QQQQC; 


01/1948 


- 12/1988 


28 


Tan-Shui (Fu-Shan, Taiwan) 


154554910 


01/1953 


- 12/1993 


29 


Pei-Nan Chi (Yen-Ping, Taiwan) 


308124037 


01/1956 


- 12/1993 


30 


Nan-Ao (Shan-Chiao, Taiwan) 


51564903 


01/1954 


- 12/1993 


31 


Cho-Shui (Tung-Tou, Taiwan) 


165580533 


01/1956 


- 12/1993 


32 


Leigh (Mount Mercer, Australia) 


134463 


10/1956 


- 12/2009 


33 


Broken (Moorngag, Australia) 


188197 


12/1973 


- 12/2009 


34 


Barwon (Inverleigh, Australia) 


290761 


04/1966 


- 12/2009 



Table 1: List of rivers (at the specified station) in the northeastern Gulf of Mexico (NEGOM) region (k = 1 — 17) and around 
the globe (k = 18 — 34), whose historical average daily discharge data is used for the analyses in this work. Data is provided 
from U.S. Geological Survey (k = 1 — 19), National River Flow Archive (k = 20 — 21), Norwegian Water Resources and 
Energy Directorate (k = 22 - 24), Royal Irrigation Department-GAME-T2 Data Center (k = 25 - 26), RosHydromet (k = 27), 
Pacific Rim Streamflow Data Set (k = 28 — 31), Department of Natural Resources and Environment- Victorian Water Resources 
(k = 32 - 34). 



ing the existence of a dynamics common to all rivers in 
the time domain, with no regards for climate, hydrologi- 
cal characteristics or geographical location. 

To better capture the universal features of the underly- 
ing dynamics we adopt a more general approach, originally 
developed for the investigation of turbulent flows at differ- 
ent time scales [TO] . It consists in analyzing a set of N new 
time series, each constructed as the difference between the 
original time series, locally detrended to remove nonsta- 
tionarity, and a copy of it delayed by n (n = 1,2, ...,N) 
days. The detrending procedure is as follows. Each time 
series is divided into sliding segments of size 2n and, for 
each segment, the trend is obtained from the best \x— th 
order polynomial P K (t) = uq + a\t + otit 2 + ... + fit- 
ting the data, where the value of (i is varied from up 



to 20. We have verified that for the time scales of inter- 
ests to our study, smaller than six years, the 20th order 
is enough to safely remove seasonal trends. As for the 
order /i of the polynomial model that best approximates 
the original time series q K (t), we choose the one that min- 
imizes the Akaike information criterion |20j . preventing 
over-fitting. The resulting time series of the difference be- 
tween the data q K (t) and the trend P K (t) is then used for 
the analysis. We indicate by Q K (t) = q K (t) — P K (t) the 
detrended time series obtained from such a procedure. It 
is worth remarking that data has not been deseasonal- 
ized before applying this procedure, because such a gen- 
eral approach is able to remove any trend, including the 
seasonal one, by construction. However, we have verified 
that pre-processing the data by the standard deseason- 



p-3 



M. Dc Domcnico V. Latora 










b 















1993 1994 1995 



1998 1999 2000 2001 



< 




Time (year) 



-30 -20 -10 10 20 30 



< 




-30 -20 -10 10 20 30 



Figure 1: (a) Historical time series corresponding to ten years 
of daily discharge for the Alabama river as measured at the 
Millers Ferry station, AL. (b) Delayed increments of the de- 
trended time series for a time delay n — 2 days. 



Figure 3: Probability density of increments, rescaled to stan- 
dard deviations a nK , for the daily discharge as measured at 
the k — th station (k = 1,2, 17) in Tab.[T] Two representa- 
tive distributions, for n = 4 days (left panel) and n = 512 days 
(right panel), are shown. 
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Figure 2: Daily discharge historical data of the k— th station 
(k = 1,2, ...,17) in Tab.[T] (a) Autocorrelation function vs. 
delay time of deseasonalized time series; (6) Power spectrum 
vs. frequency: dashed lines correspond to straight lines with 
slope -2; (c) Rank-ordered scaling exponents (P K (u) ~ u 13 ") 
corresponding to the slopes of the linear regions in Fig.^b), for 
log(i^/day _1 ) G [—2,-1] (diamonds), [—3,-2] (squares) and 
[—4, —3] (triangles). 
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Figure 4: Deformation of the PDF of standardized incre- 
ments across scales. (a) Vertically shifted PDFs of the 
Fish Creek river are shown for different time scales: n — 
8, 16, 32, 64, 128, 256, 512, 1024 and 2048 days (from top to bot- 
tom), together with the corresponding Castaing fitting curves 
(solid lines). A Gaussian distribution is also shown for refer- 
ence (dotted line) (b) Scale dependence of the non-Gaussian 
parameter vs. logn for Barwon (« = 34, upper diamonds), 
Inabon (k = 19, upper squares), Tan-Shui (k = 28, upper 
triangles), Tombigbee (k = 3, lower diamonds), Verdalsvass- 
draget (k = 24, lower squares) and Thames (k = 20, lower 
triangles). 



alizing technique has no impact on the final result. Let 
5 n Qn{t) = Q K (t + n) — Q K {t) be the time series of in- 
crements, at the time scale n, for measurements at the 
k — th station, and let (Qnn) and a\ K be the corresponding 
sample mean and variance, respectively. In the following 
we indicate by A n Q K (t) = S n Q K (t) — (Q nK ) the time se- 
ries of increments shifted to zero mean. As an example, 
in Fig.QJb) we report the time series of the standardized 
increment A n Q K (t)/a nK obtained for the Alabama river 
when we adopt a time delay of n = 2 days. 

In Fig. [3] we show the PDFs of standardized increments 
for the stations in the NEGOM region reported in Tab.[T] 
(k = 1,2, ...,17), and for two different values of time delay, 
namely n = 4 days (left panel) and n = 512 days (right 



panel). Strikingly the curves for different rivers all col- 
lapse into the same non-Gaussian distribution. Moreover, 
we have verified that for any value of the polynomial or- 
der /it, ranging from 13 up to 250, results do not change. 
Such a result was expected. In fact, the largest time scale 
we consider is n — 2048 days, corresponding to about 
six seasons: the procedure requires at most fx = 12 to 
deseasonalizc segments of length 2n, namely twelve sea- 
sons. In addition to this, we find that the same universal 
behavior holds at any time scale, from days to years. Sim- 
ilar anomalous behaviors have been observed in several 
other real-world processes of different nature, such as in 
the healthy human heart rate [2"T1|2"2"] , in stock-price fluctu- 
ations 23, 24 , in human behavioral organization and 
in seismic events |26j . 
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Figure 5: Deformation of the PDF of standardized increments 
across scales, (a) Same as in Fig.|4ja) in the case of the Tombig- 
bee river. (6) Scale dependence of the non-Gaussian parameter 
Xn vs. log n for all time series in Tab.[TJ as well as the ensemble 
average for each n. 



Models involving (eventually periodic) autoregrcssivc 
approaches, provide Gaussian PDFs of delayed increments 
and thus fail to correctly reproduce the data at any 
time scale. Conversely, heavy-tailed non-Gaussian PDFs 
emerge, for instance, from random multiplicative processes 
[THI27]- I n this case, the simplest approach is to consider 
two indipendent Gaussian distributions, F n and G n , de- 
pending on the time scale n, with zero mean and vari- 
ances respectively equal to er 2 and A 2 . If we define the 
increments A n Z{t) as 

A n Z(t) = Ut)e" n(t) , (1) 

where ^ n (t) ~ F n (0, a 2 ) and oj n (t) ~ G„(0, A 2 ), it is pos- 
sible to prove that the PDF is the Castaing function 

P°° /A Z\ 1 
p n (A n Z) = F n (-^\-G n (\na)d(\na), (2) 

originally introduced for a log-normal cascade model of 
fully developed turbulence [T5]. If delayed fluctuations 
A n Z arc standardized, it can be shown that c 2 = 
cxp(— A 2 ) and the function in Eq. ^ depends just on one 
parameter, namely A„. The standard Gaussian distribu- 
tion of delayed increments is attained when A„ approaches 
zero: for this reason A 2 has been commonly adopted as a 
key parameter in quantifying the non-Gaussian behavior 
of fluctuations [2Tll221l24If2^] . We have therefore fitted the 
PDF of delayed increments obtained from our detrended 
river flow data with the function in Eq. ([2]), to estimate 
the non-Gaussian parameter A 2 at different time scales n. 
As an example, the result of the fit is reported in Fig. 01a) 
and EJa) for the cases of the Fish Creek and Tombigbcc 
rivers, respectively. The non-Gaussian distribution of fluc- 
tuations is clear at every time scale and it is well fitted by 
a Castaing function with a value of A 2 which decreases 
with n. Wc have repeated the fitting procedure for each 
of the time series in Tab.[TJ 

Apparently, rivers exhibit a common scaling with no re- 
gards for drainage area, length, elevation, average annual 



precipitation or geographical position. Previous studies 
have already shown, for instance, the independence on the 
basin area or the geographical location of parameters char- 
acterizing the scaling behaviour in river flow data [13,28 . 
As a representative example, in Fig.@|b) is shown the de- 
pendence of the non-Gaussian parameter A 2 on the scale n 
for some rivers characterized by very different climate, hy- 
drological characteristics and location. For Barwon (Vic- 
toria, Australia), the average annual rainfall ranges from 
2000 mm, in the south of the Barwon Basin, to 700 mm 
in the northern region of the basin, whereas for Inabon 
(Ponce, Puerto Rico) the average annual rainfall reaches 
2500 mm. Although these two rivers have very different 
length (160 km and 28 km, respectively) and drainage 
area (1269 km 2 and 25 km 2 , respectively), they exhibit 
the same scaling behaviour: A 2 linearly depends on time 
scale n up to a certain horizon, where a crossover with a 
new regime appears. The same result is found for rivers 
corresponding to points in the bottom of Fig.^b), closer 
to the Gaussian limit. In this case the average annual pre- 
cipitation is rather different: for instance, 700 mm for the 
Thames (UK, Europe) and 1200-1500 mm for the Tombig- 
bee (Alabama, USA). 

In Fig.[5jb) we report, all together, the non-Gaussian 
parameters found for the 34 different rivers, and their av- 
erage, as a function of the logarithm of the time scale n. 
The scale dependence of A 2 is a common feature of all time 
scries analyzed. Our results clearly indicate that A 2 scales, 
for values of n up to 264 days, as the logarithm of n, rather 
than as a power law, while for n > 264 A 2 takes a constant 
value. In the Kolmogorov-Obukhov theory of turbulence 
[221132] the non-Gaussian parameter A 2 is expected to be 
linear in logn [31) . whereas a power law is expected in 
the case of Castaing model. Therefore, the fluctuations 
observed in river flows have a strong analogy with hydro- 
dynamic turbulence-like behavior, for temporal scales up 
to about 264 days. However, log-normal models are un- 
able to explain the nearly constant values of A 2 , observed 
for n > 264, which indicate scale invariance, and suggest 
a transition to a critical regime where the PDFs of fluc- 
tuations all perfectly collapse on the same non-Gaussian 
curve. For n > 264 days, correlations become persistent 
and strong fluctuations develop at any time scale, from the 
shortest scales above the horizon up to the largest ones. 

Our results are intriguing for different reasons. First, we 
have found a novel scaling in the time domain, common to 
all rivers, with no regards for their hydrological features or 
their location: such behaviour is peculiar of multiplicative 
cascade systems and breaks up above a certain horizon, 
where scalc-invariancc appears. Second, although the val- 
ues of the non-Gaussian parameter are not the same for 
all rivers, for a fixed time scale n, it is possible to capture 
a single dynamics, as shown in Fig.[5lb). The main differ- 
ence among rivers is not in the shape of A 2 , constrained 
by the underlying dynamics, but in the divergence from 
the limit of Gaussian fluctuations. Roughly speaking, the 
river flow dynamics appears to be always the same, and 
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rivers just differentiate because of the variance of non- 
Gaussian fluctuations they exhibit. It is worth remarking 
that multiplicative cascade models, suggested in recent 
studies [13l[T4l[T8], represent ideal candidates to explain 
our results. In fact, such models have been successfully 
adopted for modelling turbulent dynamics (see Ref. [35], 
and Rcfs. therein, for a comprehensive review of multi- 
plicative models for turbulence) and it has been recently 
shown that a single universal function, deduced from a 
generalized multiplicative model and depending just on 
two parameters, is able to reproduce the multifractal spec- 
tra of many rivers around the world |13j . The general- 
ized multiplicative model accounts for both the scaling in 
the structure function for time scales greater than several 
weeks [13] and the deformation of the PDF of standard- 
ized increment across scales up to about 264 days. The 
scale-invariance observed above the horizon cannot be ac- 
counted for when only multiplicative cascade models are 
considered [23]. Recent studies reported a crossover time 
of several weeks, above which daily runoffs are long-term 
correlated [T31fT4] . Such crossover time is about half of the 
value we have found within the present study for the tran- 
sition between two different regimes. However, the two re- 
sults can not be directly compared because they have been 
obtained from the investigation of quite different features 
of river flow, the scaling in the structure function, in one 
case, and the scaling of non-Gaussian curves in our case. 

Summing up, in this Letter we have found a novel scal- 
ing for river dynamics which relates the non-Gaussian fluc- 
tuations occurring in the flows at different time scales. 
Both the scaling behavior and the break up of the scaling 
are features common to all the rivers we have considered. 
Moreover, recent studies involving multifractal detrended 
fluctuation analysis and multifractal formalism show that 
multiplicative cascades, and their extended versions, are 
able to characterize a large number of rivers [L3l[L4lll8] . 
The universality of our results supports the appropriate- 
ness of random multiplicative models, setting up a defini- 
tive framework for modelling flow dynamics, but suggest, 
at the same time, the need for a generalization that is able 
to reproduce both well-known results in literature and the 
novel scale-invariance above the critical horizon. 

The authors thank K. Kiyono for useful discussions on 
the multiscale fluctuation analysis and the anonymous ref- 
erees for interesting comments and suggestions. 
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